/*
   \file  siagian.cc
   \brief Implementation of the siagian algorithm.

   This file defines the non-inline member functions of the
   gist::siagian class. See the header file for more info.
*/

/*
   This file is part of libgist.

   libgist is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2 of the License, or (at your
   option) any later version.

   libgist is distributed in the hope that it will be useful, but WITHOUT
   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
   FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
   for more details.

   You should have received a copy of the GNU General Public License
   along with libgist; if not, write to the Free Software Foundation,
   Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.
*/

/*
   REVISION HISTORY

   $HeadURL: http://libgist.googlecode.com/svn/branches/siagian_dev/src/siagian.cc $
   $Id: siagian.cc 142 2012-02-14 04:02:13Z christian.g.siagian@gmail.com $
*/

//------------------------------ HEADERS --------------------------------

// libgist
#include "siagian.hh"

// Boost.Thread
#include <boost/thread/locks.hpp>
#include <boost/thread/mutex.hpp>

// Boost.Timer
#include <boost/timer.hpp>

// Standard C++
#include <iostream>
#include <list>
#include <map>
#include <cstdio>
#include <vector>

// IT++ for Principal Component Analysis
#include <itpp/signal/fastica.h>

//----------------------- NAMESPACE DEFINITION --------------------------

namespace gist {

//----------------------------- CONSTANTS -------------------------------

namespace {

// constant for RGBY operation thresholding
const float RGBY_THRESH = 25.5f ;

// constants to select the proper center-surround levels 
// used by the gist features
const int delta_min = 3;
const int delta_max = 4;
const int level_min = 1;
const int level_max = 3;
const int PYR_DEPTH = (level_max + delta_max + 1);

// number of oriented pyramid levels
const  uint ORI_PYR_DEPTH = 4;

// number of reduced feature dimension
const int NUMBER_ICA_COMPONENTS = 80;

} // end of local anonymous namespace encapsulating above constants

//----------------------- GIST FEATURES ---------------------------------

// gist_features encapsulate all the gist feature values 
// extracted from the image
namespace {

// basically a vector of raw gist feature values
class gist_features {

public:

  // total number of features
  // 21 per feature maps
  // there are 12 maps in the color (rg and by) channel
  //            6 maps in the intensity channel
  //           16 maps in the orientation channel (4 angles)
  static const int SIZE = 21*(12+6+16); 

   // Create gist features 
  gist_features(std::vector<float> vals) ;

  // get the feature values
  std::vector<float> get_feature_values() ;

   // This type facilitates easy iteration over each of the SIFT
   // descriptors in the grid.
   //typedef boost::multi_array_ref<const sift_descriptor, 1> grid_data ;

   // This method returns an instance of the above type to allow clients
   // to iterate over each of the SIFT descriptors contained by this
   // grid.
   //grid_data data() const {
   //   return grid_data(m_grid.data(), boost::extents[m_grid.num_elements()]) ;
   //}

private:
  std::vector<float> its_vals;
} ;

// Constructor: create an array of specified size
gist_features::gist_features(std::vector<float> vals)
{
  its_vals.clear();
  its_vals.resize(SIZE);
  for(int i = 0; i < SIZE; i++)
    its_vals[i] = vals[i];
}

// get the feature values
std::vector<float> gist_features::get_feature_values()
{
  return its_vals;
}

// get RG and BY color opponency maps 
void getRGBY(const boost::gil::rgb8_image_t& img,
             boost::gil::gray32f_image_t& rg,
             boost::gil::gray32f_image_t& by,
             float thresh)
{
  boost::gil::rgb8_image_t::const_view_t img_view = boost::gil::const_view(img);

  //const int w = img.width() ;
  //const int h = img.height();
  //std::cout << "image size: " << w << " " << h << "\n";

  //const int wrg = rg.width() ;
  //const int hrg = rg.height();
  //const int wby = by.width() ;
  //const int hby = by.height();
  //std::cout << "image size rg: " << wrg << " " << hrg << " " << thresh <<"\n";
  //std::cout << "image size by: " << wby << " " << hby << " " << thresh <<"\n";

  //assert(w > 0 && h > 0);
  //assert(wrg == wby && hrg == hby);

  boost::gil::gray8c_step_view_t red_view   = boost::gil::nth_channel_view(img_view, 0);
  boost::gil::gray8c_step_view_t green_view = boost::gil::nth_channel_view(img_view, 1);
  boost::gil::gray8c_step_view_t blue_view  = boost::gil::nth_channel_view(img_view, 2);

  boost::gil::rgb8c_view_t::iterator img_it      = img_view.begin();
  boost::gil::rgb8c_view_t::iterator img_it_stop = img_view.end();
  
  boost::gil::gray8c_step_view_t::iterator redv_it   = red_view.begin();
  boost::gil::gray8c_step_view_t::iterator greenv_it = green_view.begin();
  boost::gil::gray8c_step_view_t::iterator bluev_it  = blue_view.begin();
  
  boost::gil::gray32f_view_t::iterator rg_it = boost::gil::view(rg).begin();
  boost::gil::gray32f_view_t::iterator by_it = boost::gil::view(by).begin();
  
  // red = [r - (g+b)/2]        [.] = clamp between 0 and 255
  // green = [g - (r+b)/2]
  // blue = [b - (r+g)/2]
  // yellow = [2*((r+g)/2 - |r-g| - b)]

  float zero = 0.0F; float thresh3 = 3.0F * thresh;
  
  while (img_it != img_it_stop)
    {
      float r = *redv_it++;
      float g = *greenv_it++;
      float b = *bluev_it++;

      // first do the luminanceNormalization:
      float lum = float(r + g + b);
      if (lum < thresh3)  // too dark... no response from anybody
        { *rg_it++ = zero; *by_it++ = zero; }
      else
        {
          // normalize chroma by luminance:
          float fac = 255.0f / lum;
          r *= fac; g *= fac; b *= fac;

          // now compute color opponencies:
          // yellow gets a factor 2 to compensate for its previous attenuation
          // by luminanceNormalize():
          float red  = r - 0.5f * (g + b), green  = g - 0.5f * (r + b);
          float blue = b - 0.5f * (r + g), yellow = -2.0f * (blue + fabs(r-g));

          if (red < 0.0f) red = 0.0f;
          else if (red > 255.0f) red = 255.0f;
          if (green < 0.0f) green = 0.0f;
          else if (green > 255.0f) green = 255.0f;
          if (blue < 0.0f) blue = 0.0f;
          else if (blue > 255.0f) blue = 255.0f;
          if (yellow < 0.0f) yellow = 0.0f;
          else if (yellow > 255.0f) yellow=255.0f;

          *rg_it++ = red - green; *by_it++ = blue - yellow;
        }
      img_it++;
    }
}

// calculate the luminance of the image
void getLuminance(const boost::gil::rgb8_image_t& img,
                  boost::gil::gray32f_image_t& lu     )
{
  boost::gil::rgb8_image_t::const_view_t img_view = boost::gil::const_view(img);

  //const int w = img.width() ;
  //const int h = img.height();
  //std::cout << "image size: " << w << " " << h << "\n";
  //assert(w > 0 && h > 0);

  boost::gil::gray8c_step_view_t red_view   = boost::gil::nth_channel_view(img_view, 0);
  boost::gil::gray8c_step_view_t green_view = boost::gil::nth_channel_view(img_view, 1);
  boost::gil::gray8c_step_view_t blue_view  = boost::gil::nth_channel_view(img_view, 2);

  //boost::gil::rgb8c_view_t::iterator img_it      = img_view.begin();
  //boost::gil::rgb8c_view_t::iterator img_it_stop = img_view.end();
  
  boost::gil::gray8c_step_view_t::iterator redv_it   = red_view.begin();
  boost::gil::gray8c_step_view_t::iterator greenv_it = green_view.begin();
  boost::gil::gray8c_step_view_t::iterator bluev_it  = blue_view.begin();
  boost::gil::gray8c_step_view_t::iterator redv_stop = red_view.end();
  
  boost::gil::gray32f_view_t::iterator lu_it = boost::gil::view(lu).begin();

  while (redv_it != redv_stop)
    {
      // float r = *redv_it++;
      // float g = *greenv_it++;
      // float b = *bluev_it++;
      // float lum = (r + g + b)/3.0f;
      // *lu_it++ = lum;
      *lu_it++ = (*redv_it++ + *greenv_it++ + *bluev_it++)/3.0f;
    }
}

// Anderson's separable kernel: 1/16 * [1 4 6 4 1]
boost::gil::gray32f_image_t
lowPass5yDecY(boost::gil::gray32f_image_t& source)
{
  // assume factor = 2
  const uint factor = 2;

  //assert(factor > 0);

  // if factor == 1 then we don't decimate:
  //if (factor == 1) return lowPass5y(src);

  // now factor is guaranteed to be >= 2
  const int w = source.width(), h = source.height();
  int h2 = h / factor;
  if (h2 == 0) h2 = 1; // do not go smaller than 1 pixel

  //if (h < 2) return source; // nothing to smooth

  boost::gil::gray32f_image_t result(w,h2);

  boost::gil::gray32f_view_t::iterator sptr = boost::gil::view(source).begin();
  boost::gil::gray32f_view_t::iterator dptr = boost::gil::view(result).begin();

  // ########## vertical pass  (even though we scan horiz for speedup)
  const int w2 = w * 2, w3 = w * 3, w4 = w * 4; // speedup

  if (h == 2) //////////////////////////////////////////////////
    {
      // topmost points  ( [ (6^) 4 ] / 10 )^T
      for (int i = 0; i < w; i ++)
        {
          //*dptr++ = sptr[0] * (6.0F / 10.0F) + sptr[w] * (4.0F / 10.0F);
          float v0 = boost::gil::at_c<0>(sptr[0]);
          float v1 = boost::gil::at_c<0>(sptr[w]);
          *dptr++ = v0 * (6.0F / 10.0F) + v1 * (4.0F / 10.0F);
          sptr++;
        }
      sptr -= w;  // go back to top-left
    }
  else if (h == 3) //////////////////////////////////////////////////
    {
      // topmost points  ( [ (6^) 4 1 ] / 11 )^T
      for (int i = 0; i < w; i ++)
        {
          // *dptr++ = sptr[ 0] * (6.0F / 11.0F) +
          //   sptr[ w] * (4.0F / 11.0F) +
          //   sptr[w2] * (1.0F / 11.0F);
          float v0 = boost::gil::at_c<0>(sptr[0 ]);
          float v1 = boost::gil::at_c<0>(sptr[w ]);
          float v2 = boost::gil::at_c<0>(sptr[w2]);
          *dptr++ = v0 * (6.0F / 11.0F) + v1 * (4.0F / 11.0F) + v2 * (1.0F / 11.0F);
          sptr++;
        }
      sptr -= w;  // go back to top-left

    // since factor >= 2 we don't need any other point.
  }
  else  ///////////////////////////////// general case for height >= 4
    {
      int i1 = 0, i2 = 0;
      const int skip = (factor - 1) * w;

      // topmost points  ( [ (6^) 4 1 ] / 11 )^T
      for (int i = 0; i < w; i ++)
        {
          //*dptr++ = sptr[ 0] * (6.0F / 11.0F) +
          //  sptr[ w] * (4.0F / 11.0F) +
          //  sptr[w2] * (1.0F / 11.0F);
          float v0 = boost::gil::at_c<0>(sptr[0 ]);
          float v1 = boost::gil::at_c<0>(sptr[w ]);
          float v2 = boost::gil::at_c<0>(sptr[w2]);
          *dptr++ = v0 * (6.0F / 11.0F) + v1 * (4.0F / 11.0F) + v2 * (1.0F / 11.0F);
          sptr++;
        }
      sptr -= w;  // go back to top-left
      ++i2;
      i1 += factor;

      // second point skipped since factor >= 2
      sptr += (skip - w);

      // rest of the column except last 2 points ( [ 1^ 4 (6) 4 1 ] / 16 )T
      while((i1 < (h-2)) && (i2 < h2))
        {
          for (int i = 0; i < w; i ++)
            {
              // *dptr++ = (sptr[ 0] + sptr[w4]) * (1.0F / 16.0F) +
              //   (sptr[ w] + sptr[w3]) * (4.0F / 16.0F) +
              //   sptr[w2]  * (6.0F / 16.0F);
              float v0 = boost::gil::at_c<0>(sptr[0 ]);
              float v1 = boost::gil::at_c<0>(sptr[w ]);
              float v2 = boost::gil::at_c<0>(sptr[w2]);
              float v3 = boost::gil::at_c<0>(sptr[w3]);
              float v4 = boost::gil::at_c<0>(sptr[w4]);
              *dptr++ = (v0 + v4) * (1.0F / 16.0F) + 
                (v1 + v3) * (4.0F / 16.0F) + 
                v2  * (6.0F / 16.0F);
              sptr++;
            }
          sptr += skip;
          i1 += factor;
          ++ i2;
        }

      // need special case for second to last point?
      if ((i2 < h2) && (i1 == (h-2)))
        {
          sptr = boost::gil::view(source).end() - w4;
          // before last points ( [ 1^ 4 (6) 4 ] / 15 )T
          for (int i = 0; i < w; i ++)
            {
              // *dptr++ = sptr[ 0] * (1.0F / 15.0F) +
              //   (sptr[ w] + sptr[w3]) * (4.0F / 15.0F) +
              //   sptr[w2] * (6.0F / 15.0F);
              float v0 = boost::gil::at_c<0>(sptr[0 ]);
              float v1 = boost::gil::at_c<0>(sptr[w ]);
              float v2 = boost::gil::at_c<0>(sptr[w2]);
              float v3 = boost::gil::at_c<0>(sptr[w3]);
              *dptr++ = v0 * (1.0F / 15.0F) +      
                (v1 + v3) * (4.0F / 15.0F) + 
                v2 * (6.0F / 15.0F);
              sptr++;
            }
          i1 += factor;
          ++i2;
        }

      // need special case for last point?
      if ((i2 < h2) && (i1 == (h-1)))
        {
          sptr = boost::gil::view(source).end() - w3;
          // last points ( [ 1^ 4 (6) ] / 11 )T
          for (int i = 0; i < w; i ++)
            {
              // *dptr++ = sptr[ 0] * (1.0F / 11.0F) +
              //   sptr[ w] * (4.0F / 11.0F) +
              //   sptr[w2] * (6.0F / 11.0F);
              float v0 = boost::gil::at_c<0>(sptr[0 ]);
              float v1 = boost::gil::at_c<0>(sptr[w ]);
              float v2 = boost::gil::at_c<0>(sptr[w2]);
              *dptr++ = v0 * (1.0F / 11.0F) +
                v1 * (4.0F / 11.0F) +
                v2 * (6.0F / 11.0F);
              sptr++;
            }
        }
    }

  return result;
}

// Anderson's separable kernel: 1/16 * [1 4 6 4 1]
boost::gil::gray32f_image_t
lowPass5xDecX(boost::gil::gray32f_image_t& source)
{
  // assume factor = 2
  const uint factor = 2;

  assert(factor > 0);

  // if factor == 1 then we don't decimate:
  //if (factor == 1) return lowPass5x(src);

  // now factor is guaranteed to be >= 2
  const int w = source.width(), h = source.height();
  int w2 = w / factor;
  if (w2 == 0) w2 = 1; // do not go smaller than 1 pixel

  //if (w < 2) return source; // nothing to smooth

  boost::gil::gray32f_image_t result(w2,h);

  boost::gil::gray32f_view_t::iterator sptr = boost::gil::view(source).begin();
  boost::gil::gray32f_view_t::iterator dptr = boost::gil::view(result).begin();


  if (w == 2) //////////////////////////////////////////////////
    for (int j = 0; j < h; j ++)
      {
        // leftmost point  [ (6^) 4 ] / 10
        //*dptr++ = sptr[0] * (6.0F / 10.0F) + sptr[1] * (4.0F / 10.0F);
        float v0 = boost::gil::at_c<0>(sptr[0]);
        float v1 = boost::gil::at_c<0>(sptr[1]);
        *dptr++ = v0 * (6.0F / 10.0F) + v1 * (4.0F / 10.0F);

        sptr += 2;  // sptr back to same position as dptr
      }
  else if (w == 3) //////////////////////////////////////////////////
    for (int j = 0; j < h; j ++)
      {
        // need the left most point in any case
        // leftmost point  [ (6^) 4 1 ] / 11
        // *dptr++ = sptr[0] * (6.0F / 11.0F) +
        //   sptr[1] * (4.0F / 11.0F) +
        //   sptr[2] * (1.0F / 11.0F);
        float v0 = boost::gil::at_c<0>(sptr[0 ]);
        float v1 = boost::gil::at_c<0>(sptr[1 ]);
        float v2 = boost::gil::at_c<0>(sptr[2 ]);
        *dptr++ = v0 * (6.0F / 11.0F) +
          v1 * (4.0F / 11.0F) +
          v2 * (1.0F / 11.0F);

        // since factor >= 2 we don't need any other point.

        sptr += 3;  // sptr back to same position as dptr
      }
  else  ////////////////////////////// general case for width() >= 4
    // *** unfolded version (all particular cases treated) for
    // max speed.  do not use access overloading since can
    // access the C array directly.  the bet is: float
    // computations are faster than int (true for most float
    // copros), and int2float convert is fast -> do not have to
    // divide by the coeff sum since we directly use float
    // coeffs.
    // notations: in () is the position of dest ptr, and ^ is src ptr
    // ########## horizontal pass
    for (int j = 0; j < h; j ++)
      {
        int i1 = 0, i2 = 0;
        boost::gil::gray32f_view_t::iterator sptr2 = sptr;

        // leftmost point  [ (6^) 4 1 ] / 11
        // *dptr++ = sptr2[0] * (6.0F / 11.0F) +
        //   sptr2[1] * (4.0F / 11.0F) +
        //   sptr2[2] * (1.0F / 11.0F);
        float v20 = boost::gil::at_c<0>(sptr2[0 ]);
        float v21 = boost::gil::at_c<0>(sptr2[1 ]);
        float v22 = boost::gil::at_c<0>(sptr2[2 ]);
        *dptr++ = v20 * (6.0F / 11.0F) +
          v21 * (4.0F / 11.0F) +
          v22 * (1.0F / 11.0F);
        ++i2;
        i1 += factor;

        // skip second point since factor >= 2:
        sptr2 += (factor-2);
        
        // rest of the line except last 2 points  [ 1^ 4 (6) 4 1 ] / 16.0
        while ((i1 < (w-2)) && (i2 < w2))
          {
            // *dptr++ = (sptr2[0] + sptr2[4]) * (1.0F / 16.0F) +
            //   (sptr2[1] + sptr2[3]) * (4.0F / 16.0F) +
            //   sptr2[2]  * (6.0F / 16.0F);
            float v20 = boost::gil::at_c<0>(sptr2[0 ]);
            float v21 = boost::gil::at_c<0>(sptr2[1 ]);
            float v22 = boost::gil::at_c<0>(sptr2[2 ]);
            float v23 = boost::gil::at_c<0>(sptr2[3 ]);
            float v24 = boost::gil::at_c<0>(sptr2[4 ]);
            *dptr++ = (v20 + v24) * (1.0F / 16.0F) +
              (v21 + v23) * (4.0F / 16.0F) +
              v22  * (6.0F / 16.0F);
            i1 += factor; sptr2 += factor;
            ++i2;
          }

        // need special case for second to last point?
        if ((i2 < w2) && (i1 == (w-2)))
          {
            sptr2 = sptr + w - 4;
            // before last point [ 1^ 4 (6) 4 ] / 15
            // *dptr++ = sptr2[0] * (1.0F / 15.0F) +
            //   (sptr2[1] + sptr2[3]) * (4.0F / 15.0F) +
            //   sptr2[2] * (6.0F / 15.0F);
            float v20 = boost::gil::at_c<0>(sptr2[0 ]);
            float v21 = boost::gil::at_c<0>(sptr2[1 ]);
            float v22 = boost::gil::at_c<0>(sptr2[2 ]);
            float v23 = boost::gil::at_c<0>(sptr2[3 ]);
            *dptr++ = v20 * (1.0F / 15.0F) +
              (v21 + v23) * (4.0F / 15.0F) +
              v22 * (6.0F / 15.0F);
            i1 += factor;
            ++i2;
          }

        // need special case for last point?
        if ((i2 < w2) && (i1 == (w-1)))
          {
            sptr2 = sptr + w - 3;
            // last point [ 1^ 4 (6) ] / 11
            // *dptr++ = sptr2[0] * (1.0F / 11.0F) +
            //   sptr2[1] * (4.0F / 11.0F) +
            //   sptr2[2] * (6.0F / 11.0F);
            float v20 = boost::gil::at_c<0>(sptr2[0 ]);
            float v21 = boost::gil::at_c<0>(sptr2[1 ]);
            float v22 = boost::gil::at_c<0>(sptr2[2 ]);
            *dptr++ = v20 * (1.0F / 11.0F) +
              v21 * (4.0F / 11.0F) +
              v22 * (6.0F / 11.0F);
            ++i2;
          }
        sptr += w;
      }

  return result;
}

// helper printing function
void print(boost::gil::gray32f_image_t& a)
{
  const int w = a.width(), h = a.height();
  printf("[%d %d]\n",w,h);

  boost::gil::gray32f_view_t av = boost::gil::view(a);
  boost::gil::gray32f_view_t::xy_locator av_loc = av.xy_at(0,0);
  for(uint i = 0; i < 3; i++)
    {
      for(uint j = 0; j < 3; j++)
        {
          float aval = boost::gil::at_c<0>(av_loc(i,j));
          printf("[%d %d]: %f\n", i, j, aval);
        }
    }
}

// seperable convolution for smoothing maps
boost::gil::gray32f_image_t sepFilterX9(boost::gil::gray32f_image_t source)
{
  // assume filter size = 9
  int hfs = 9;

  // assume boundary policy: CONV_BOUNDARY_CLEAN

  float hFilt[9] = {  1.0F / 256.0F,  8.0F / 256.0F, 28.0F / 256.0F,
                     56.0F / 256.0F, 70.0F / 256.0F, 56.0F / 256.0F,
                     28.0F / 256.0F,  8.0F / 256.0F,  1.0F / 256.0F };
  
  // assume filter_size = 9
  boost::gil::gray32f_image_t result(source.dimensions());

  const int w = source.width(), h = source.height();

  //if (hfs == 0) return source;  // no filter

  boost::gil::gray32f_view_t::iterator sptr = boost::gil::view(source).begin();
  boost::gil::gray32f_view_t::iterator dptr = boost::gil::view(result).begin();

  const int hfs2 = (hfs - 1) / 2;

  // flip the filter to accelerate convolution:
  float *hFilter = new float[hfs]; float sumh = 0.0f;
  for (int x = 0; x < hfs; x ++)
    { sumh += hFilt[x]; hFilter[hfs-1-x] = hFilt[x]; }

  // *** horizontal pass ***
  if (w < hfs)  // special function for very small images
    {
      for (int j = 0; j < h; j ++)
        for (int i = 0; i < w; i ++)
          {
            float val = 0.0F; float sum = 0.0F;
            for (int k = 0; k < hfs; k ++)
              if (i + k - hfs2 >= 0 && i + k - hfs2 < w)
                {
                  val += boost::gil::at_c<0>(sptr[k - hfs2]) * hFilter[k];
                  sum += hFilter[k];
                }
            *dptr++ = val * sumh / sum; sptr ++;
          }
    }
  else  // function for reasonably large images
    for (int jj = 0; jj < h; jj ++)
      {
        // leftmost points:
        for (int j = hfs2; j > 0; --j)
          {
            float val = 0.0f; float sum = 0.0F;
            for (int k = j; k < hfs; ++k)
              {
                val += boost::gil::at_c<0>(sptr[k - j]) * hFilter[k];
                sum += hFilter[k];
              }
            *dptr++ = val * sumh / sum;
          }

        // bulk (far from the borders):
        for (int i = 0; i < w - hfs + 1; i ++)
          {
            float val = 0.0f;
            for (int k = 0; k < hfs; k ++) 
              val += boost::gil::at_c<0>(sptr[k]) * hFilter[k];
            *dptr++ = val; sptr++;
          }

        // rightmost points:
        for (int j = hfs - 1; j > hfs2; --j)
          {
            float val = 0.0f; float sum = 0.0F;
            for (int k = 0; k < j; ++k)
              {
                val += boost::gil::at_c<0>(sptr[k]) * hFilter[k];
                sum += hFilter[k];
              }
            *dptr++ = val * sumh / sum;
            sptr++;
          }

        sptr += hfs2;  // sptr back to same as dptr (start of next line)
      }
  //delete [] hFilter;
  return result;
}

// separable convolution for smoothing maps
boost::gil::gray32f_image_t sepFilterY9(boost::gil::gray32f_image_t source)
{
  // assume filter size = 9
  int vfs = 9;

  // assume boundary policy: CONV_BOUNDARY_CLEAN

   float vFilt[9] = {  1.0F / 256.0F, 8.0F / 256.0F, 28.0F / 256.0F,
                      56.0F / 256.0F, 70.0F / 256.0F, 56.0F / 256.0F,
                      28.0F / 256.0F, 8.0F / 256.0F, 1.0F / 256.0F };

  boost::gil::gray32f_image_t result(source.dimensions());

  const int w = source.width(), h = source.height();

  //if (vfs == 0) return source;  // no filter

  boost::gil::gray32f_view_t::iterator sptr = boost::gil::view(source).begin();
  boost::gil::gray32f_view_t::iterator dptr = boost::gil::view(result).begin();

  const int vfs2 = (vfs - 1) / 2;

  // flip the filter to accelerate convolution:
  float *vFilter = new float[vfs]; float sumv = 0.0f;
  for (int x = 0; x < vfs; x ++)
    { sumv += vFilt[x]; vFilter[vfs-1-x] = vFilt[x]; }

  int ww[vfs];
  for (int i = 0; i < vfs; i ++) ww[i] = w * i; // speedup precompute

  if (h < vfs)  // special function for very small images
    {
      for (int j = 0; j < h; j ++)
        for (int i = 0; i < w; i ++)
          {
            float val = 0.0f; float sum = 0.0;
            for (int k = 0; k < vfs; k ++)
              if (j + k - vfs2 >= 0 && j + k - vfs2 < h)
                {
                  val += boost::gil::at_c<0>(sptr[w * (k - vfs2)]) * vFilter[k];
                  sum += vFilter[k];
                }
            *dptr++ = val * sumv / sum; sptr ++;
          }
    }
  else  // function for reasonably large images
    {
      // top points:
      for (int j = vfs2; j > 0; --j)
        {
          for (int i = 0; i < w; ++i)  // scan all points of given horiz
            {
              float val = 0.0f; float sum = 0.0;
              for (int k = j; k < vfs; ++k)
                {
                  val += boost::gil::at_c<0>(sptr[ww[k - j]]) * vFilter[k];
                  sum += vFilter[k];
                }
              *dptr++ = val * sumv / sum;
              sptr++;
            }
          sptr -= w;   // back to top-left corner
        }

      // bulk (far from edges):
      for (int j = 0; j < h - vfs + 1; j ++)
        for (int i = 0; i < w; i ++)
          {
            float val = 0.0f;
            for (int k = 0; k < vfs; k ++) 
              val += boost::gil::at_c<0>(sptr[ww[k]]) * vFilter[k];
            *dptr++ = val; sptr ++;
          }

      // bottommost points:
      for (int j = vfs - 1; j > vfs2; --j)
        for (int i = 0; i < w; ++i)
          {
            float val = 0.0f; float sum = 0.0;
            for (int k = 0; k < j; ++k)
              {
                val += boost::gil::at_c<0>(sptr[ww[k]]) * vFilter[k];
                sum += vFilter[k];
              }
            *dptr++ = val * sumv / sum;
            sptr ++;
          }
    }
  return result;
}

// low pass the image in the X direction
boost::gil::gray32f_image_t lowpassX(boost::gil::gray32f_image_t source)
{
  const int w = source.width(), h = source.height();
  //if (w < 2) return source; // nothing to smooth

  boost::gil::gray32f_image_t result(source.dimensions());

  boost::gil::gray32f_view_t::iterator sptr = boost::gil::view(source).begin();
  boost::gil::gray32f_view_t::iterator dptr = boost::gil::view(result).begin();

  if (w < 9)  // use inefficient implementation for small images
    {
      return sepFilterX9(source);
    }

  // boundary conditions: truncated filter
  for (int j = 0; j < h; j ++)
    {
      float v0 = boost::gil::at_c<0>(sptr[0 ]);
      float v1 = boost::gil::at_c<0>(sptr[1 ]);
      float v2 = boost::gil::at_c<0>(sptr[2 ]);
      float v3 = boost::gil::at_c<0>(sptr[3 ]);
      float v4 = boost::gil::at_c<0>(sptr[4 ]);
      float v5 = boost::gil::at_c<0>(sptr[5 ]);
      float v6 = boost::gil::at_c<0>(sptr[6 ]);
      float v7 = boost::gil::at_c<0>(sptr[7 ]);
      float v8 = boost::gil::at_c<0>(sptr[8 ]);

      // leftmost points
      // *dptr++ = sptr[0] * (70.0F / 163.0F) +
      //   sptr[1] * (56.0F / 163.0F) +
      //   sptr[2] * (28.0F / 163.0F) +
      //   sptr[3] * ( 8.0F / 163.0F) +
      //   sptr[4] * ( 1.0F / 163.0F);
      *dptr++ = v0 * (70.0F / 163.0F) +
        v1 * (56.0F / 163.0F) +
        v2 * (28.0F / 163.0F) +
        v3 * ( 8.0F / 163.0F) +
        v4 * ( 1.0F / 163.0F);

      // *dptr++ = (sptr[0] + sptr[2]) * (56.0F / 219.0F) +
      //   sptr[1] * (70.0F / 219.0F) +
      //   sptr[3] * (28.0F / 219.0F) +
      //   sptr[4] * ( 8.0F / 219.0F) +
      //   sptr[5] * ( 1.0F / 219.0F);
      *dptr++ = (v0 + v2) * (56.0F / 219.0F) +
        v1 * (70.0F / 219.0F) +
        v3 * (28.0F / 219.0F) +
        v4 * ( 8.0F / 219.0F) +
        v5 * ( 1.0F / 219.0F);

      // *dptr++ = (sptr[0] + sptr[4]) * (28.0F / 247.0F) +
      //   (sptr[1] + sptr[3]) * (56.0F / 247.0F) +
      //   sptr[2] * (70.0F / 247.0F) +
      //   sptr[5] * ( 8.0F / 247.0F) +
      //   sptr[6] * ( 1.0F / 247.0F);
      *dptr++ = (v0 + v4) * (28.0F / 247.0F) +
        (v1 + v3) * (56.0F / 247.0F) +
        v2 * (70.0F / 247.0F) +
        v5 * ( 8.0F / 247.0F) +
        v6 * ( 1.0F / 247.0F);

      // *dptr++ = (sptr[0] + sptr[6]) * ( 8.0F / 255.0F) +
      //   (sptr[1] + sptr[5]) * (28.0F / 255.0F) +
      //   (sptr[2] + sptr[4]) * (56.0F / 255.0F) +
      //   sptr[3] * (70.0F / 255.0F) +
      //   sptr[7] * ( 1.0F / 255.0F);
      *dptr++ = (v0 + v6) * ( 8.0F / 255.0F) +
        (v1 + v5) * (28.0F / 255.0F) +
        (v2 + v4) * (56.0F / 255.0F) +
        v3 * (70.0F / 255.0F) +
        v7 * ( 1.0F / 255.0F);

      // far from the borders
      for (int i = 0; i < w - 8; i ++)
        {
          v0 = boost::gil::at_c<0>(sptr[0 ]);
          v1 = boost::gil::at_c<0>(sptr[1 ]);
          v2 = boost::gil::at_c<0>(sptr[2 ]);
          v3 = boost::gil::at_c<0>(sptr[3 ]);
          v4 = boost::gil::at_c<0>(sptr[4 ]);
          v5 = boost::gil::at_c<0>(sptr[5 ]);
          v6 = boost::gil::at_c<0>(sptr[6 ]);
          v7 = boost::gil::at_c<0>(sptr[7 ]);
          v8 = boost::gil::at_c<0>(sptr[8 ]);

          // *dptr++ = (sptr[0] + sptr[8]) * ( 1.0F / 256.0F) +
          //   (sptr[1] + sptr[7]) * ( 8.0F / 256.0F) +
          //   (sptr[2] + sptr[6]) * (28.0F / 256.0F) +
          //   (sptr[3] + sptr[5]) * (56.0F / 256.0F) +
          //   sptr[4] * (70.0F / 256.0F);
          *dptr++ = (v0 + v8) * ( 1.0F / 256.0F) +
            (v1 + v7) * ( 8.0F / 256.0F) +
            (v2 + v6) * (28.0F / 256.0F) +
            (v3 + v5) * (56.0F / 256.0F) +
            v4 * (70.0F / 256.0F);
          sptr ++;
        }

      v0 = boost::gil::at_c<0>(sptr[0 ]);
      v1 = boost::gil::at_c<0>(sptr[1 ]);
      v2 = boost::gil::at_c<0>(sptr[2 ]);
      v3 = boost::gil::at_c<0>(sptr[3 ]);
      v4 = boost::gil::at_c<0>(sptr[4 ]);
      v5 = boost::gil::at_c<0>(sptr[5 ]);
      v6 = boost::gil::at_c<0>(sptr[6 ]);
      v7 = boost::gil::at_c<0>(sptr[7 ]);
      v8 = boost::gil::at_c<0>(sptr[8 ]);

      // rightmost points
      // *dptr++ = sptr[0] * ( 1.0F / 255.0F) +
      //   (sptr[1] + sptr[7]) * ( 8.0F / 255.0F) +
      //   (sptr[2] + sptr[6]) * (28.0F / 255.0F) +
      //   (sptr[3] + sptr[5]) * (56.0F / 255.0F) +
      //   sptr[4] * (70.0F / 255.0F);
      *dptr++ = v0 * ( 1.0F / 255.0F) +
        (v1 + v7) * ( 8.0F / 255.0F) +
        (v2 + v6) * (28.0F / 255.0F) +
        (v3 + v5) * (56.0F / 255.0F) +
        v4 * (70.0F / 255.0F);
      sptr ++;

      v0 = boost::gil::at_c<0>(sptr[0 ]);
      v1 = boost::gil::at_c<0>(sptr[1 ]);
      v2 = boost::gil::at_c<0>(sptr[2 ]);
      v3 = boost::gil::at_c<0>(sptr[3 ]);
      v4 = boost::gil::at_c<0>(sptr[4 ]);
      v5 = boost::gil::at_c<0>(sptr[5 ]);
      v6 = boost::gil::at_c<0>(sptr[6 ]);
      v7 = boost::gil::at_c<0>(sptr[7 ]);
      v8 = boost::gil::at_c<0>(sptr[8 ]);

      // *dptr++ = sptr[0] * ( 1.0F / 247.0F) +
      //   sptr[1] * ( 8.0F / 247.0F) +
      //   (sptr[2] + sptr[6]) * (28.0F / 247.0F) +
      //   (sptr[3] + sptr[5]) * (56.0F / 247.0F) +
      //   sptr[4] * (70.0F / 247.0F);
      *dptr++ = v0 * ( 1.0F / 247.0F) +
        v1 * ( 8.0F / 247.0F) +
        (v2 + v6) * (28.0F / 247.0F) +
        (v3 + v5) * (56.0F / 247.0F) +
        v4 * (70.0F / 247.0F);
      sptr ++;

      v0 = boost::gil::at_c<0>(sptr[0 ]);
      v1 = boost::gil::at_c<0>(sptr[1 ]);
      v2 = boost::gil::at_c<0>(sptr[2 ]);
      v3 = boost::gil::at_c<0>(sptr[3 ]);
      v4 = boost::gil::at_c<0>(sptr[4 ]);
      v5 = boost::gil::at_c<0>(sptr[5 ]);
      v6 = boost::gil::at_c<0>(sptr[6 ]);
      v7 = boost::gil::at_c<0>(sptr[7 ]);
      v8 = boost::gil::at_c<0>(sptr[8 ]);

      // *dptr++ = sptr[0] * ( 1.0F / 219.0F) +
      //   sptr[1] * ( 8.0F / 219.0F) +
      //   sptr[2] * (28.0F / 219.0F) +
      //   (sptr[3] + sptr[5]) * (56.0F / 219.0F) +
      //   sptr[4] * (70.0F / 219.0F);
      *dptr++ = v0 * ( 1.0F / 219.0F) +
        v1 * ( 8.0F / 219.0F) +
        v2 * (28.0F / 219.0F) +
        (v3 + v5) * (56.0F / 219.0F) +
        v4 * (70.0F / 219.0F);

      sptr ++;

      v0 = boost::gil::at_c<0>(sptr[0 ]);
      v1 = boost::gil::at_c<0>(sptr[1 ]);
      v2 = boost::gil::at_c<0>(sptr[2 ]);
      v3 = boost::gil::at_c<0>(sptr[3 ]);
      v4 = boost::gil::at_c<0>(sptr[4 ]);
      v5 = boost::gil::at_c<0>(sptr[5 ]);
      v6 = boost::gil::at_c<0>(sptr[6 ]);
      v7 = boost::gil::at_c<0>(sptr[7 ]);
      v8 = boost::gil::at_c<0>(sptr[8 ]);

      // *dptr++ = sptr[0] * ( 1.0F / 163.0F) +
      //   sptr[1] * ( 8.0F / 163.0F) +
      //   sptr[2] * (28.0F / 163.0F) +
      //   sptr[3] * (56.0F / 163.0F) +
      //   sptr[4] * (70.0F / 163.0F);
      *dptr++ = v0 * ( 1.0F / 163.0F) +
        v1 * ( 8.0F / 163.0F) +
        v2 * (28.0F / 163.0F) +
        v3 * (56.0F / 163.0F) +
        v4 * (70.0F / 163.0F);

      sptr += 5;  // sptr back to same as dptr (start of next line)
    }

  return result;
}

// low pass the image
boost::gil::gray32f_image_t lowpassY(boost::gil::gray32f_image_t source)
{
  // assume filter_size = 9
  const int w = source.width(), h = source.height();
  //if (h < 2) return source; // nothing to smooth

  boost::gil::gray32f_image_t result(source.dimensions());

  boost::gil::gray32f_view_t::iterator sptr = boost::gil::view(source).begin();
  boost::gil::gray32f_view_t::iterator dptr = boost::gil::view(result).begin();

  if (h < 9)  // use inefficient implementation for small images
    {
      return sepFilterY9(source);
    }

  // *** vertical pass ***
  const int w2 = w + w, w3 = w2 + w, w4 = w3 + w, w5 = w4 + w, w6 = w5 + w,
    w7 = w6 + w,  w8 = w7 + w;  // index computation speedup
  for (int i = 0; i < w; i ++)
    {
      float v0 = boost::gil::at_c<0>(sptr[0 ]);
      float v1 = boost::gil::at_c<0>(sptr[w ]);
      float v2 = boost::gil::at_c<0>(sptr[w2]);
      float v3 = boost::gil::at_c<0>(sptr[w3]);
      float v4 = boost::gil::at_c<0>(sptr[w4]);

      // *dptr++ = sptr[ 0] * (70.0F / 163.0F) +
      //   sptr[ w] * (56.0F / 163.0F) +
      //   sptr[w2] * (28.0F / 163.0F) +
      //   sptr[w3] * ( 8.0F / 163.0F) +
      //   sptr[w4] * ( 1.0F / 163.0F);
      *dptr++ = v0 * (70.0F / 163.0F) +
        v1 * (56.0F / 163.0F) +
        v2 * (28.0F / 163.0F) +
        v3 * ( 8.0F / 163.0F) +
        v4 * ( 1.0F / 163.0F);
      sptr ++;
    }
  sptr -= w; // back to top-left
  for (int i = 0; i < w; i ++)
    {
      float v0 = boost::gil::at_c<0>(sptr[0 ]);
      float v1 = boost::gil::at_c<0>(sptr[w ]);
      float v2 = boost::gil::at_c<0>(sptr[w2]);
      float v3 = boost::gil::at_c<0>(sptr[w3]);
      float v4 = boost::gil::at_c<0>(sptr[w4]);
      float v5 = boost::gil::at_c<0>(sptr[w5]);

      // *dptr++ = (sptr[ 0] + sptr[w2]) * (56.0F / 219.0F) +
      //   sptr[ w] * (70.0F / 219.0F) +
      //   sptr[w3] * (28.0F / 219.0F) +
      //   sptr[w4] * ( 8.0F / 219.0F) +
      //   sptr[w5] * ( 1.0F / 219.0F);
      *dptr++ = (v0 + v2) * (56.0F / 219.0F) +
        v1 * (70.0F / 219.0F) +
        v3 * (28.0F / 219.0F) +
        v4 * ( 8.0F / 219.0F) +
        v5 * ( 1.0F / 219.0F);

      sptr ++;
    }
  sptr -= w; // back to top-left
  for (int i = 0; i < w; i ++)
    {
      float v0 = boost::gil::at_c<0>(sptr[0 ]);
      float v1 = boost::gil::at_c<0>(sptr[w ]);
      float v2 = boost::gil::at_c<0>(sptr[w2]);
      float v3 = boost::gil::at_c<0>(sptr[w3]);
      float v4 = boost::gil::at_c<0>(sptr[w4]);
      float v5 = boost::gil::at_c<0>(sptr[w5]);
      float v6 = boost::gil::at_c<0>(sptr[w6]);

      // *dptr++ = (sptr[ 0] + sptr[w4]) * (28.0F / 247.0F) +
      //   (sptr[ w] + sptr[w3]) * (56.0F / 247.0F) +
      //   sptr[w2] * (70.0F / 247.0F) +
      //   sptr[w5] * ( 8.0F / 247.0F) +
      //   sptr[w6] * ( 1.0F / 247.0F);
      *dptr++ = (v0 + v4) * (28.0F / 247.0F) +
        (v1 + v3) * (56.0F / 247.0F) +
        v2 * (70.0F / 247.0F) +
        v5 * ( 8.0F / 247.0F) +
        v6 * ( 1.0F / 247.0F);
      sptr ++;
    }
  sptr -= w; // back to top-left
  for (int i = 0; i < w; i ++)
    {
      float v0 = boost::gil::at_c<0>(sptr[0 ]);
      float v1 = boost::gil::at_c<0>(sptr[w ]);
      float v2 = boost::gil::at_c<0>(sptr[w2]);
      float v3 = boost::gil::at_c<0>(sptr[w3]);
      float v4 = boost::gil::at_c<0>(sptr[w4]);
      float v5 = boost::gil::at_c<0>(sptr[w5]);
      float v6 = boost::gil::at_c<0>(sptr[w6]);
      float v7 = boost::gil::at_c<0>(sptr[w7]);

      // *dptr++ = (sptr[ 0] + sptr[w6]) * ( 8.0F / 255.0F) +
      //   (sptr[ w] + sptr[w5]) * (28.0F / 255.0F) +
      //   (sptr[w2] + sptr[w4]) * (56.0F / 255.0F) +
      //   sptr[w3] * (70.0F / 255.0F) +
      //   sptr[w7] * ( 1.0F / 255.0F);
      *dptr++ = (v0 + v6) * ( 8.0F / 255.0F) +
        (v1 + v5) * (28.0F / 255.0F) +
        (v2 + v4) * (56.0F / 255.0F) +
        v3 * (70.0F / 255.0F) +
        v7 * ( 1.0F / 255.0F);
      sptr ++;
    }
  sptr -= w;   // back to top-left
  for (int j = 0; j < h - 8; j ++)
    for (int i = 0; i < w; i ++)
      {
        float v0 = boost::gil::at_c<0>(sptr[0 ]);
        float v1 = boost::gil::at_c<0>(sptr[w ]);
        float v2 = boost::gil::at_c<0>(sptr[w2]);
        float v3 = boost::gil::at_c<0>(sptr[w3]);
        float v4 = boost::gil::at_c<0>(sptr[w4]);
        float v5 = boost::gil::at_c<0>(sptr[w5]);
        float v6 = boost::gil::at_c<0>(sptr[w6]);
        float v7 = boost::gil::at_c<0>(sptr[w7]);
        float v8 = boost::gil::at_c<0>(sptr[w8]);

        // *dptr++ = (sptr[ 0] + sptr[w8]) * ( 1.0F / 256.0F) +
        //   (sptr[ w] + sptr[w7]) * ( 8.0F / 256.0F) +
        //   (sptr[w2] + sptr[w6]) * (28.0F / 256.0F) +
        //   (sptr[w3] + sptr[w5]) * (56.0F / 256.0F) +
        //   sptr[w4]  * (70.0F / 256.0F);
        *dptr++ = (v0 + v8) * ( 1.0F / 256.0F) +
          (v1 + v7) * ( 8.0F / 256.0F) +
          (v2 + v6) * (28.0F / 256.0F) +
          (v3 + v5) * (56.0F / 256.0F) +
          v4  * (70.0F / 256.0F);
        sptr ++;
      }
  for (int i = 0; i < w; i ++)
    {
      float v0 = boost::gil::at_c<0>(sptr[0 ]);
      float v1 = boost::gil::at_c<0>(sptr[w ]);
      float v2 = boost::gil::at_c<0>(sptr[w2]);
      float v3 = boost::gil::at_c<0>(sptr[w3]);
      float v4 = boost::gil::at_c<0>(sptr[w4]);
      float v5 = boost::gil::at_c<0>(sptr[w5]);
      float v6 = boost::gil::at_c<0>(sptr[w6]);
      float v7 = boost::gil::at_c<0>(sptr[w7]);

      // *dptr++ = sptr[ 0] * ( 1.0F / 255.0F) +
      //   (sptr[ w] + sptr[w7]) * ( 8.0F / 255.0F) +
      //   (sptr[w2] + sptr[w6]) * (28.0F / 255.0F) +
      //   (sptr[w3] + sptr[w5]) * (56.0F / 255.0F) +
      //   sptr[w4] * (70.0F / 255.0F);
      *dptr++ = v0 * ( 1.0F / 255.0F) +
        (v1 + v7) * ( 8.0F / 255.0F) +
        (v2 + v6) * (28.0F / 255.0F) +
        (v3 + v5) * (56.0F / 255.0F) +
        v4 * (70.0F / 255.0F);
      sptr ++;
    }
  for (int i = 0; i < w; i ++)
    {
      float v0 = boost::gil::at_c<0>(sptr[0 ]);
      float v1 = boost::gil::at_c<0>(sptr[w ]);
      float v2 = boost::gil::at_c<0>(sptr[w2]);
      float v3 = boost::gil::at_c<0>(sptr[w3]);
      float v4 = boost::gil::at_c<0>(sptr[w4]);
      float v5 = boost::gil::at_c<0>(sptr[w5]);
      float v6 = boost::gil::at_c<0>(sptr[w6]);

      // *dptr++ = sptr[ 0] * ( 1.0F / 247.0F) +
      //   sptr[ w] * ( 8.0F / 247.0F) +
      //   (sptr[w2] + sptr[w6]) * (28.0F / 247.0F) +
      //   (sptr[w3] + sptr[w5]) * (56.0F / 247.0F) +
      //   sptr[w4] * (70.0F / 247.0F);
      *dptr++ = v0 * ( 1.0F / 247.0F) +
        v1 * ( 8.0F / 247.0F) +
        (v2 + v6) * (28.0F / 247.0F) +
        (v3 + v5) * (56.0F / 247.0F) +
        v4 * (70.0F / 247.0F);
      sptr ++;
    }
  for (int i = 0; i < w; i ++)
    {
      float v0 = boost::gil::at_c<0>(sptr[0 ]);
      float v1 = boost::gil::at_c<0>(sptr[w ]);
      float v2 = boost::gil::at_c<0>(sptr[w2]);
      float v3 = boost::gil::at_c<0>(sptr[w3]);
      float v4 = boost::gil::at_c<0>(sptr[w4]);
      float v5 = boost::gil::at_c<0>(sptr[w5]);

      // *dptr++ = sptr[ 0] * ( 1.0F / 219.0F) +
      //   sptr[ w] * ( 8.0F / 219.0F) +
      //   sptr[w2] * (28.0F / 219.0F) +
      //   (sptr[w3] + sptr[w5]) * (56.0F / 219.0F) +
      //   sptr[w4] * (70.0F / 219.0F);
      *dptr++ = v0 * ( 1.0F / 219.0F) +
        v1 * ( 8.0F / 219.0F) +
        v2 * (28.0F / 219.0F) +
        (v3 + v5) * (56.0F / 219.0F) +
        v4 * (70.0F / 219.0F);
      sptr ++;
    }
  for (int i = 0; i < w; i ++)
    {
      float v0 = boost::gil::at_c<0>(sptr[0 ]);
      float v1 = boost::gil::at_c<0>(sptr[w ]);
      float v2 = boost::gil::at_c<0>(sptr[w2]);
      float v3 = boost::gil::at_c<0>(sptr[w3]);
      float v4 = boost::gil::at_c<0>(sptr[w4]);

      // *dptr++ = sptr[ 0] * ( 1.0F / 163.0F) +
      //   sptr[ w] * ( 8.0F / 163.0F) +
      //   sptr[w2] * (28.0F / 163.0F) +
      //   sptr[w3] * (56.0F / 163.0F) +
      //   sptr[w4] * (70.0F / 163.0F);
      *dptr++ = v0 * ( 1.0F / 163.0F) +
        v1 * ( 8.0F / 163.0F) +
        v2 * (28.0F / 163.0F) +
        v3 * (56.0F / 163.0F) +
        v4 * (70.0F / 163.0F);
      sptr ++;
    }

  return result;
}

// low pass the image
boost::gil::gray32f_image_t lowPass(uint filter_size, boost::gil::gray32f_image_t image)
{
  //assume filter_size = 9
  filter_size += 0;

  boost::gil::gray32f_image_t result(image);

  result = lowpassX(result); 
  result = lowpassY(result); 

  return result;
}

// FIXXX quick subtraction function, should be done more elequently
boost::gil::gray32f_image_t
subtract(boost::gil::gray32f_image_t& a,
         boost::gil::gray32f_image_t& b)
{
  boost::gil::gray32f_image_t r(a.dimensions());

  boost::gil::gray32f_view_t::iterator a_it = boost::gil::view(a).begin();
  boost::gil::gray32f_view_t::iterator b_it = boost::gil::view(b).begin();
  boost::gil::gray32f_view_t::iterator r_it = boost::gil::view(r).begin();
  boost::gil::gray32f_view_t::iterator a_it_stop = 
    boost::gil::view(a).end();

  while (a_it != a_it_stop)
    {
      // FIXXX there has got to be a better way
      //*r_it++ = *a_it++ - *b_it++;
      *r_it++ = 
        boost::gil::at_c<0>(*a_it++) - boost::gil::at_c<0>(*b_it++);      
    }
  return r;
}

// decrementing the image, right now only to half size
boost::gil::gray32f_image_t
decXY(boost::gil::gray32f_image_t& src)
{
  // assume decrement factor of 2
  //const int yfactor = yfactor_raw >= 0 ? yfactor_raw : xfactor;
  const int xfactor = 2;
  const int yfactor = 2;

  // do not go smaller than 1x1:
  //if (src.width() <= 1 && src.height() <= 1) return src;

  if (src.width() == 1)  // only thinout vertic
    {
      int h2 = src.height() / yfactor;
      boost::gil::gray32f_image_t result(1,h2);

      boost::gil::gray32f_view_t::iterator a_it = 
        boost::gil::view(src).begin();
      boost::gil::gray32f_view_t::iterator r_it = 
        boost::gil::view(result).begin();
      boost::gil::gray32f_view_t::iterator r_it_stop = 
        boost::gil::view(result).end();
      
      uint j = 0;
      while (r_it != r_it_stop)
        { *r_it++ = a_it[j*yfactor]; j++; }

      return result;
    }
  if (src.height() == 1)
    {
      int w2 = src.width() / xfactor;
      boost::gil::gray32f_image_t result(w2,1);

      boost::gil::gray32f_view_t::iterator a_it = 
        boost::gil::view(src).begin();
      boost::gil::gray32f_view_t::iterator r_it = 
        boost::gil::view(result).begin();
      boost::gil::gray32f_view_t::iterator r_it_stop = 
        boost::gil::view(result).end();
      
      uint i = 0;
      while (r_it != r_it_stop)
        { *r_it++ = a_it[i*xfactor]; i++; }

      return result;
    }

  const int w2 = src.width()  / xfactor;
  const int h2 = src.height() / yfactor;

  boost::gil::gray32f_image_t result(w2,h2);

  boost::gil::gray32f_view_t::iterator sptr = 
    boost::gil::view(src).begin();
  boost::gil::gray32f_view_t::iterator dptr = 
    boost::gil::view(result).begin();

  const int skip = src.width() % xfactor + src.width() * (yfactor - 1);

  for (int j = 0; j < h2; ++j)
    {
      for (int i = 0; i < w2; ++i)
        {
          *dptr++ = *sptr;   // copy one pixel
          sptr += xfactor;   // skip some pixels
        }
      sptr += skip;          // skip to start of next line
    }

  return result;
}

// get laplacian pyramid
std::vector<boost::gil::gray32f_image_t> 
buildPyrLaplacian(boost::gil::gray32f_image_t& image, uint depth, 
                  uint filter_size)
{
  // assume first level = 0

  std::vector<boost::gil::gray32f_image_t> result(depth); 
  boost::gil::gray32f_image_t lpf = lowPass(filter_size, image);
  
  //if (0 >= firstlevel)
  result[0] = subtract(image,lpf);

  for (uint lev = 1; lev < depth; lev++)
    {
      boost::gil::gray32f_image_t dec = decXY(lpf);
      lpf = lowPass(filter_size, dec);

      //if (lev >= firstlevel)
      result[lev] = subtract(dec,lpf);
    }

  return result;
}

// filter image in particular orientation
boost::gil::gray32f_image_t orientedFilter
(boost::gil::gray32f_image_t src, 
 float spatial_freq, float theta, float intens)
{
  boost::gil::gray32f_image_t result(src.dimensions());

  // assume usetab = false

  const float k = spatial_freq;
  const float intensity = intens;

  double kx = double(k) * cos((theta + 90.0) * M_PI / 180.0);
  double ky = double(k) * sin((theta + 90.0) * M_PI / 180.0);

  boost::gil::gray32f_image_t re(src.dimensions());
  boost::gil::gray32f_image_t im(src.dimensions());

  boost::gil::gray32f_view_t::iterator sptr = 
    boost::gil::view(src).begin();
  boost::gil::gray32f_view_t::iterator reptr = 
    boost::gil::view(re).begin();
  boost::gil::gray32f_view_t::iterator imptr = 
    boost::gil::view(im).begin();

  // (x,y) = (0,0) at center of image:
  int w2l = src.width() / 2 , w2r = src.width() - w2l;
  int h2l = src.height() / 2, h2r = src.height() - h2l;

  for (int j = -h2l; j < h2r; ++j)
    for (int i = -w2l; i < w2r; ++i)
      {
        const double arg = kx * i + ky * j;
        const float val = boost::gil::at_c<0>(*sptr++) * intensity;

        // Otherwise we resort to explicit separate sin() and
        // cos() calls.
        double sinarg = sin(arg);
        double cosarg = cos(arg);

        *reptr++ = float(val * cosarg);
        *imptr++ = float(val * sinarg);

      }

  re = lowpassX(re);
  re = lowpassY(re);
   
  im = lowpassX(im);
  im = lowpassY(im);

  // quad energy
  boost::gil::gray32f_view_t::iterator dptr = 
    boost::gil::view(result).begin();
  boost::gil::gray32f_view_t::iterator stop = 
    boost::gil::view(result).end();
  boost::gil::gray32f_view_t::iterator s1ptr = 
    boost::gil::view(re).begin();
  boost::gil::gray32f_view_t::iterator s2ptr = 
    boost::gil::view(im).begin();

  while (dptr != stop)
    {
      const float s1( boost::gil::at_c<0>(*s1ptr++) );
      const float s2( boost::gil::at_c<0>(*s2ptr++) );
      *dptr++ = float(sqrt(s1 * s1 + s2 * s2));
    }

  return result;
}

// border attenuation
void inplaceAttenuateBorders(boost::gil::gray32f_image_t& a, 
                             float attenuation_width)
{
  boost::gil::gray32f_image_t result(a.dimensions());

  int size = attenuation_width;  

  int w = a.width();
  int h = a.height();

  if (size * 2 > w) size = w / 2;
  if (size * 2 > h) size = h / 2;
  if (size < 1) return;  // forget it

  float increment = 1.0 / (float)(size + 1);
  // top lines:
  float coeff = increment;
  boost::gil::gray32f_view_t::iterator aptr = boost::gil::view(a).begin();

  for (int y = 0; y < size; y ++)
    {
      for (int x = 0; x < w; x ++)
        {
          *aptr = boost::gil::at_c<0>(*aptr) * coeff;
          ++aptr;
        }
      coeff += increment;
    }
  // normal lines: start again from beginning to attenuate corners twice:
  aptr =  boost::gil::view(a).begin();
  for (int y = 0; y < h; y ++)
    {
      coeff = increment;
      for (int x = 0; x < size; x ++)
        {
          *(aptr + w - 1 - x * 2) =
            boost::gil::at_c<0>(*(aptr + w - 1 - x * 2)) * coeff;

          *aptr = boost::gil::at_c<0>(*aptr) * coeff;
          ++aptr;
          coeff += increment;
        }
      aptr += w - size;
    }
  // bottom lines
  aptr =  boost::gil::view(a).begin() + (h - size) * w;
  coeff = increment * (float)size;
  for (int y = h - size; y < h; y ++)
    {
      for (int x = 0; x < w; ++x)
        {
          *aptr = boost::gil::at_c<0>(*aptr) * coeff;
          ++aptr;
        }
      coeff -= increment;
    }
}

// build oriented pyramid from laplacian
std::vector<boost::gil::gray32f_image_t> 
buildPyrOrientedFromLaplacian
(std::vector<boost::gil::gray32f_image_t>& laplacian, 
 float theta, float intens)
{
  // assume first level = 0
  // assume filter size = 9

  int attenuation_width = 5;
  float spatial_freq = 2.6;

  uint depth = laplacian.size();

  // compute oriented filter from laplacian:
  std::vector<boost::gil::gray32f_image_t> result(depth); 

  for (uint lev = 0; lev < result.size(); ++lev)
    {
      // if (laplacian[lev].initialized())
      // {
      //result[lev] = laplacian[lev];
      result[lev] = orientedFilter(laplacian[lev], spatial_freq,
                                   theta, intens);

      // attenuate borders that are overestimated due to filter trunctation:
      inplaceAttenuateBorders(result[lev], attenuation_width);

      // }
    }
  return result;
}

// get gaussian pyramid
std::vector<boost::gil::gray32f_image_t> 
getGaussianPyramid(boost::gil::gray32f_image_t& ima, uint depth)
{
  // assume first level = 0
  // assume filter size = 5

  std::vector<boost::gil::gray32f_image_t> result(depth); 
  boost::gil::gray32f_image_t a(ima);
  result[0] = a;  
  
  for (uint lev = 1; lev < depth; lev++)
    {
      a = lowPass5xDecX(a); a = lowPass5yDecY(a);
      //a = lowPass5yDecY(lowPass5xDecX(a,2),2); can't do it in 1 line
      result[lev] = a;
    }

  return result;
}

// get oriented pyramid
std::vector<boost::gil::gray32f_image_t> 
getOrientedPyramid(boost::gil::gray32f_image_t& ima, uint depth, float theta)
{
  // assume first level = 0
  // assume filter size = 9
  const uint filter_size = 9;
  // assume intens = 10.0f
  const float intens = 10.0f;

  std::vector<boost::gil::gray32f_image_t> laplacian =
    buildPyrLaplacian(ima, depth, filter_size);

  std::vector<boost::gil::gray32f_image_t>  orifilt =
    buildPyrOrientedFromLaplacian(laplacian, theta, intens);

  return orifilt;
}

// get center surround maps
std::vector<boost::gil::gray32f_image_t> 
getCenterSurround(std::vector<boost::gil::gray32f_image_t> pyr)
{
  std::vector<boost::gil::gray32f_image_t> cs; 

  // assume absol = true

  for (int lev = level_min; lev <= level_max; lev ++)
    for (int delta = delta_min; delta <= delta_max; delta ++)
      {
        boost::gil::gray32f_image_t center   = pyr[lev      ];
        boost::gil::gray32f_image_t surround = pyr[lev+delta];

        const int lw = center.width()  , lh = center.height();
        const int sw = surround.width(), sh = surround.height();

        int scalex = lw / sw, remx = lw - 1 - (lw % sw);
        int scaley = lh / sh, remy = lh - 1 - (lh % sh);

        // result has the size of the larger image:
        boost::gil::gray32f_image_t result(center.dimensions());
        
        // scan large image and subtract corresponding pixel from small image:
        int ci = 0, cj = 0;
        boost::gil::gray32f_view_t::iterator lptr = 
          boost::gil::view(center).begin();
        boost::gil::gray32f_view_t::iterator sptr = 
          boost::gil::view(surround).begin();
        boost::gil::gray32f_view_t::iterator dptr = 
          boost::gil::view(result).begin();

        for (int j = 0; j < lh; ++j)
          {
            for (int i = 0; i < lw; ++i)
              {
                if ( boost::gil::at_c<0>(*lptr) >  boost::gil::at_c<0>(*sptr))
                  *dptr++ =  boost::gil::at_c<0>(*lptr++) - boost::gil::at_c<0>(*sptr);
                else
                  *dptr++ =  boost::gil::at_c<0>(*sptr)   - boost::gil::at_c<0>(*lptr++);

                if ((++ci) == scalex && i != remx) { ci = 0; ++sptr; }
              }
            if (ci) { ci = 0; ++sptr; }  // in case the reduction is not round
            if ((++cj) == scaley && j != remy) cj = 0; else sptr -= sw;
          }

        // attenuate borders:
        int dim = result.width(); if(result.height() > dim) dim = result.height();
         inplaceAttenuateBorders(result, dim/20);
         cs.push_back(result);
      }

  return cs;
}

// get the raw gist features using subsum algorithm
std::vector<float> 
getSubSum(boost::gil::gray32f_image_t img)
{
  std::vector<float> res(21); 
  
  int w = img.width();  int h = img.height();
  std::vector<float> tempRes(21);
  for (int i = 0; i < 21; i++) tempRes[i] = 0.0;
  std::vector<int> counts(16);
  for (int i = 0 ; i < 16; i++) counts[i] = 0;

  boost::gil::gray32f_view_t::iterator itr = boost::gil::view(img).begin();
  for (int y = 0; y < h; ++y)
    {
      int suby = (4*y)/h;
      for (int x = 0; x < w; ++x)
        {
          int subx = (4*x)/w;
          int subpos = 4*suby + subx;
          tempRes[subpos+5] += boost::gil::at_c<0>(*itr);

          ++(counts[subpos]);
          ++itr;
        }

    }

  int order[] = { 5,6,9,10, 7,8,11,12, 13,14,17,18, 15,16,19,20 };
  for (int i = 0 ; i < 16; ++i)
    if (counts[i] > 0)
      res[i+5] =  tempRes[order[i]] / (counts[order[i] - 5]+ 0.0);

  float tre1 = tempRes[5] + tempRes[6] + tempRes[9] + tempRes[10];
  int ct1 = counts[0] + counts[1] + counts[4] + counts[5];

  float tre2 = tempRes[7] + tempRes[8] + tempRes[11] + tempRes[12];
  int ct2 = counts[2] + counts[3] + counts[6] + counts[7];

  float tre3 = tempRes[13] + tempRes[14] + tempRes[17] + tempRes[18];
  int ct3 = counts[8] + counts[9] + counts[12] + counts[13];

  float tre4 = tempRes[15] + tempRes[16] + tempRes[19] + tempRes[20];
  int ct4 = counts[10] + counts[11] + counts[14] + counts[15];

  res[1] = tre1/ct1;
  res[2] = tre2/ct2;
  res[3] = tre3/ct3;
  res[4] = tre4/ct4;

  res[0] = (tre1 + tre2 + tre3 + tre4)/(ct1 + ct2 + ct3 + ct4);

  // printf("lev1   : %14.4f \n", res[0]);
  // printf("lev2   : %14.4f, %14.4f\n", res[1], res[2]);
  // printf("       : %14.4f, %14.4f\n", res[3], res[4]);
  // printf("lev3   : %14.4f, %14.4f, %14.4f, %14.4f\n",
  //        res[ 5], res[ 6], res[ 9], res[10]);
  // printf("       : %14.4f, %14.4f, %14.4f, %14.4f\n",
  //        res[ 7], res[ 8], res[11], res[12]);
  // printf("       : %14.4f, %14.4f, %14.4f, %14.4f\n",
  //        res[13], res[14], res[17], res[18]);
  // printf("       : %14.4f, %14.4f, %14.4f, %14.4f\n",
  //        res[15], res[16], res[19], res[20]);

  return res;
}

// given an input image, compute its gist features
gist_features compute_gist_features(const boost::gil::rgb8_image_t& img)
{
  //const int w = img.width() ;
  //const int h = img.height();
  //std::cout << "image size: " << w << " " << h << "\n";

  // create RG, BY, and luminance images
  boost::gil::gray32f_image_t rg(img.dimensions());
  boost::gil::gray32f_image_t by(img.dimensions());
  boost::gil::gray32f_image_t lu(img.dimensions());
  getRGBY(img, rg, by, RGBY_THRESH);
  getLuminance(img, lu);

  // create gaussian and oriented pyramid builder
  std::vector<boost::gil::gray32f_image_t> rg_pyr = 
    getGaussianPyramid(rg, PYR_DEPTH);
  std::vector<boost::gil::gray32f_image_t> by_pyr = 
    getGaussianPyramid(by, PYR_DEPTH);
  std::vector<boost::gil::gray32f_image_t> in_pyr = 
    getGaussianPyramid(lu, PYR_DEPTH);

  std::vector<boost::gil::gray32f_image_t> o0_pyr   = 
    getOrientedPyramid(lu, ORI_PYR_DEPTH,   0.0f);  
  std::vector<boost::gil::gray32f_image_t> o45_pyr  = 
    getOrientedPyramid(lu, ORI_PYR_DEPTH,  45.0f);  
  std::vector<boost::gil::gray32f_image_t> o90_pyr  = 
    getOrientedPyramid(lu, ORI_PYR_DEPTH,  90.0f);  
  std::vector<boost::gil::gray32f_image_t> o135_pyr = 
    getOrientedPyramid(lu, ORI_PYR_DEPTH, 135.0f);  


  // compute center surround for rg, by, intensity sub-channel
  std::vector<boost::gil::gray32f_image_t> rg_cs = getCenterSurround(rg_pyr); 
  std::vector<boost::gil::gray32f_image_t> by_cs = getCenterSurround(by_pyr); 
  std::vector<boost::gil::gray32f_image_t> in_cs = getCenterSurround(in_pyr); 

  // get the raw gist features
  std::vector<float> gist_vals;  
  for(uint i = 0; i < rg_cs.size(); i++)
    {
      std::vector<float> vals = getSubSum(rg_cs[i]); 
      for(uint j = 0; j < vals.size(); j++)
        gist_vals.push_back(vals[j]);
    }

  for(uint i = 0; i < by_cs.size(); i++)
    {      
      std::vector<float> vals = getSubSum(by_cs[i]); 
      for(uint j = 0; j < vals.size(); j++)
        gist_vals.push_back(vals[j]);
    }

  for(uint i = 0; i < rg_cs.size(); i++)
    {
      std::vector<float> vals = getSubSum(in_cs[i]); 
      for(uint j = 0; j < vals.size(); j++)
        gist_vals.push_back(vals[j]);
    }

  for(uint i = 0; i < ORI_PYR_DEPTH; i++)
    {
      std::vector<float> vals = getSubSum(o0_pyr[i]); 
      for(uint j = 0; j < vals.size(); j++)
        gist_vals.push_back(vals[j]);
    }

  for(uint i = 0; i < ORI_PYR_DEPTH; i++)
    {
      std::vector<float> vals = getSubSum(o45_pyr[i]); 
      for(uint j = 0; j < vals.size(); j++)
        gist_vals.push_back(vals[j]);
    }

  for(uint i = 0; i < ORI_PYR_DEPTH; i++)
    {
      std::vector<float> vals = getSubSum(o90_pyr[i]); 
      for(uint j = 0; j < vals.size(); j++)
        gist_vals.push_back(vals[j]);
    }

  for(uint i = 0; i < ORI_PYR_DEPTH; i++)
    {
      std::vector<float> vals = getSubSum(o135_pyr[i]); 
      for(uint j = 0; j < vals.size(); j++)
        gist_vals.push_back(vals[j]);
    }

  // store to gist features
  gist_features g(gist_vals);

  return g ;
}

} // end of local anonymous namespace encapsulating above helpers

//----------------------- GIST FEATURES LISTS ---------------------------

namespace {

// This data structure holds all the gist features for a particular image
// category. It is used in training mode.
typedef std::list<gist_features> gist_features_list ;

// Quick helper to create a new list of gist features containing
// the grid passed in to the function.
inline gist_features_list new_list(const gist_features& g)
{
   gist_features_list list ;
   list.push_back(g) ;
   return list ;
}

} // end of local anonymous namespace encapsulating above helpers


//---------------------- SIAGIAN IMPLEMENTATION ------------------------

// The siagian class hides its entire implementation behind an opaque
// object. This section defines this implementation class.
namespace {

// In training mode, we will accumulate the gist features for the
// training images in a map indexed by the category labels.
typedef std::map<std::string, gist_features_list> gist_features_map ;

// This class actually implements the Siagian algorithm. The interface
// class defined in the header file, for the most part, simply forwards
// its API to this one.
class siagian_impl {

   // This data structure holds gist features for the input
   // images arranged according to their categories.
   gist_features_map gf_map ;

   // As we allow multithreaded clients, we will have to synchronize
   // access to the above members.
   mutable boost::mutex m_mutex ;

public:
   // Constructor
   siagian_impl() ;

   // Add a new gist features given its category
   void add(const std::string& cat, const gist_features&) ;

  gist_features_map get_gist_features_map();

} ;

// Constructor
siagian_impl::siagian_impl()
{}

// Lock map of gist features and add new features given its category
void siagian_impl::add(const std::string& cat, const gist_features& gfeat)
{
  boost::lock_guard<boost::mutex> lock(m_mutex) ;
  gist_features_map::iterator m = gf_map.find(cat) ;
  if (m == gf_map.end())
    gf_map.insert(std::make_pair(cat, new_list(gfeat))) ;
  else
    m->second.push_back(gfeat) ;
}

gist_features_map siagian_impl::get_gist_features_map()
{
  // FIXXX probably not the most efficient!!!
  return gf_map;
}

} // end of local anonymous namespace encapsulating above implementation details

//------------------------ SIAGIAN INTERFACE ---------------------------

// Initialization
siagian::siagian()
   : p_impl(new siagian_impl)
{}

// Training
void
siagian::
data(const std::string& cat, const boost::gil::rgb8_image_t& img)
{
  // compute and store the gist features
  reinterpret_cast<siagian_impl*>(p_impl)->add(cat, compute_gist_features(img)) ;
}

void siagian::train()
{
  boost::timer t ;   
  siagian_impl* pimpl = reinterpret_cast<siagian_impl*>(p_impl) ;
  double s = t.elapsed();
  std::cout << "data taking " << s << " seconds\n" ;

   // FIXXX: Lazebnik does local helper copy "class copy"
  
  //std::cout << "please enter a character\n"; 
  //char str[256]; std::cin >> str;


  std::cout << "\nSIAGIAN: PCA and ICA dimension reduction\n";

  gist_features_map gf_map = pimpl->get_gist_features_map() ;  
  uint num_labels = gf_map.size();  
  uint num_examples = 0;
  gist_features_map::iterator gptr = gf_map.begin();  
  for (uint l = 0; l < num_labels; ++l) 
    {      
      num_examples += (*gptr++).second.size();
      //std::cout << l << " " << num_examples << "\n"; 
    }

  itpp::mat  signal(gist_features::SIZE,num_examples);
  
  gptr = gf_map.begin(); uint curr_gf = 0;
  for (uint l = 0; l < num_labels; ++l) 
    {      
      std::list<gist_features>  gf_list = (*gptr++).second;
      uint size = gf_list.size();
      std::list<gist_features>::iterator lptr = gf_list.begin();
      for(uint i = 0; i < size; i++)
        {
          std::vector<float> vals = (*lptr++).get_feature_values();
          for (int c = 0; c < gist_features::SIZE; ++c) 
            {
              signal(c,curr_gf) = vals[c];
            }
          curr_gf++;
        }
    }
  double s2 = t.elapsed();
  std::cout << "moving data " << s2 << " seconds\n" ;

  std::cout << "mat: " << signal.rows() << " " << signal.cols() << " " 
            << NUMBER_ICA_COMPONENTS << "\n";  

  itpp::Fast_ICA fastICA(signal);
  fastICA.set_nrof_independent_components(NUMBER_ICA_COMPONENTS);
  // fastICA.set_non_linearity(FICA_NONLIN_TANH); // or FICA_NONLIN_POW3 
  fastICA.set_approach( FICA_APPROACH_DEFL ); // Much faster
  fastICA.separate();

  // the separating matrix
  itpp::mat W = fastICA.get_separating_matrix();

  std::cout << "W: " << W.rows() << " " << W.cols() << "\n";

  // or just the separated signals
  //itpp::mat ica_signal = fastICA.get_independent_components();
  //std::cout << "ica_signal: " << ica_signal.rows() << " " << ica_signal.cols() << "\n";

  uint n_ic = fastICA.get_nrof_independent_components();
  std::cout << "num components: " << n_ic << "\n";

  // save it

  // to run it: ica_signal = W*signal; 

  double s3 = t.elapsed();
  std::cout << "PCA/ICA takes: " << s3 << " seconds\n" ;

  // train 
  std::cout << "\nSIAGIAN: train SVM place classifier FIXXX: TODO\n";



}

// Clean-up
siagian::~siagian()
{
   delete reinterpret_cast<siagian_impl*>(p_impl) ;
}

//-----------------------------------------------------------------------

} // end of namespace encapsulating this file's definitions

/* So things look consistent in everyone's emacs... */
/* Local Variables: */
/* indent-tabs-mode: nil */
/* End: */
